Find task description and steps here
Find and download the source code from here
The model (metarange) predicts abundance of the species for all the cells in the landscape (in this case two states in Germany - Bavaria and Baden-Württemberg).
all_data <- list.files(recursive = T) %>% #get a list of the paths to datasets files
map_dfr(function(path) {
#retrieve any characters from the beginning of the path until /
foldername <- str_extract(path, "^.+(?=/)")
fread(path) %>% #import the files
mutate(folder = foldername) # add a new collumn with the folder name
})
invisible(gc())
all_data <- all_data %>%
separate(folder, c("date", "time", "species", "scenario"), sep = "_", remove = FALSE) %>% #takes a VERY long time ^(more than 40 minutes)
dplyr::select(-c("folder", "date", "time"))
invisible(gc())
#check for negative values in abundance
summary(all_data)
#rows with negative abundances
bb <- as.data.frame(which(all_data < 0 , arr.ind=TRUE))
bb <- bb[bb$col == 4,]
#values <- as.vector(bb$row)
aa <- all_data[values,]
#write tsv file with all species in all scenarios
write_tsv(all_data, "full_dataset.tsv")
#Import the complete dataset
all_data <- fread("./full_dataset.tsv", integer64 = "numeric")
nvisible(gc())Occupancy - cells with at least one individual
The following variables are all calculated in relation to time step 25, which corresponds to the end of the burn-in period that allows the species to stabilize (the model has reached it stationary region).
Abundance change - \(abundance_{tn} - abundance_{t25}\)
Occupancy change - \(ocupancy_{tn} - occupancy_{t25}\)
Habitat change - \(habitat\ suitability_{tn} - habitat\ suitability_{t25}\)
#################################
# STEP 2 - CREATE NEW VARIABLES #
#################################
# occupancy
all_data <- all_data %>%
mutate(occupancy = ifelse(all_data$abundance > 0, 1, 0))
invisible(gc())
# abundance change
all_data <- all_data %>%
group_by(species, scenario) %>%
mutate(abund_change <- abundance - abundance[t == 25]) %>%
ungroup()
invisible(gc())
# occupancy change
all_data <- all_data %>%
group_by(species, scenario) %>%
mutate(occup_change <- occupancy - occupancy[t == 25]) %>%
ungroup()
invisible(gc())
# habitat suitability change
all_data <- all_data %>%
group_by(species, scenario) %>%
mutate(habitat_change <- habitat - habitat[t == 25]) %>%
ungroup()
invisible(gc())Population-level = value for the entire population (e.g. average of cell values for the entire population in a given time step)
For the following plots the burn in period is not included in the graphical representation, since it corresponds to the preliminary steps of the model.
##########################################
# STEP 3 - CREATE POPULATION LEVEL PLOTS #
##########################################
# calculate mean values
data_all_mean_stdev<-all_data[, c(1,4:10)]
# mean and standard deviation for ABUNDANCE, REPRODUCTION, CARRYING CAPACITY, HABITAT SUITABILITY and MORTALITY variables
mean_stDEV <- data_all_mean_stdev %>%
dplyr::group_by(species, scenario, t) %>%
dplyr::summarise(across(c(abundance, reproduction, carry, habitat, bevmort), .fns = list(mean = mean, sd = sd), na.rm = TRUE), .groups = 'drop') %>%
dplyr::mutate(across(where(is.numeric), round, 3)) %>%
dplyr::filter(t>25) #remove burn-in
# write table into .csv file
write.csv(mean_stDEV, "mean_values_31Aug.csv", row.names=FALSE)
invisible(gc())abund_plot## Don't know how to automatically pick scale for object of type <integer64>.
## Defaulting to continuous.
datatable(abund_quant, caption = "Mean abundance values for 2100 (t=115) and 2011 (t=25) and quantification of the abundance change between 2100 and 2011")reprod_plothabitat_plotcarry_cap_plotmort_plotAs of 30th August 2024, the abundance mismatch was calculated only for landscape cells considered suitable (meaning it was calculated only for cells where the habitat suitability in time step 115 increased in relation the habitat suitability in time step 25). Additionally only time step 115 was plotted.**
#################################
## ABUNDANCE MISMATCH BOXPLOTS ##
#################################
abund_mismatch_df <- all_data[which(all_data$habitat_change >0), ]
invisible(gc())
abund_mismatch_df <- abund_mismatch_df[abund_mismatch_df$t==115,]
invisible(gc())
#summary(abund_mismatch_df)
# abundance mismatch
abund_mismatch_df <- abund_mismatch_df %>%
group_by(species, scenario) %>%
mutate(abund_mismatch <- carry - abundance) %>%
ungroup() %>%
rename(abund_mismatch = 15)
invisible(gc())abund_mismatch_plot_GWF_withoutliers## Don't know how to automatically pick scale for object of type <integer64>.
## Defaulting to continuous.
abund_mismatch_plot_GWF## Don't know how to automatically pick scale for object of type <integer64>.
## Defaulting to continuous.
##################################################
# STEP 4 - CREATE CELL LEVEL PLOTS (per species) #
##################################################
# abundance change proportion
all_data <- all_data %>%
group_by(species, scenario) %>%
mutate(abund_change_prop <- (abundance - abundance[t == 25])/abundance[t == 25]) %>%
rename(abund_change_prop = 15) %>%
ungroup()
invisible(gc())
# subset data for the last time step (t = 90)
abund_t115<-all_data[all_data$t == 115, c(2,3,9,10,4,15)]
#get the names of each scenario
scenarios <- unique(abund_t115$scenario)
# Define the RdBu palette with a midpoint at 0
palette <- brewer.pal(11, "RdBu")These plots were created with for loops to allow for an image with a map per species for each scenario. In total there are 7 images with 9 maps each.
###
Plot for Scenario SSP1-RCP2.6
###
Plot for Scenario SSP1
###
Plot for Scenario RCP8.5
###
Plot for Scenario SSP5-RCP8.5
###
Plot for Scenario SSP5
###
Plot for Scenario control
datatable(abund_chnage_prop_quant, caption = " Quantification of the proportion of abundance chnage between 2100 and 2011")An extra plot was created representing Spatial hotspots of change in the landscape across all 9 species (MPCAS = Mean proportion of change across species).
spatial_hotspots_mapcorrelation_plot## Warning: Removed 185535 rows containing non-finite outside the scale range
## (`stat_smooth()`).